home *** CD-ROM | disk | FTP | other *** search
/ Super PC 31 / Super PC 31 (Shareware).iso / spc / inter / speakf / fuente / lpc / lpc.c next >
Encoding:
C/C++ Source or Header  |  1995-11-07  |  6.3 KB  |  276 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "lpc.h"
  4. #define FAR _far
  5. #include "../ulaw.h"
  6. #include <stdlib.h>
  7. #include <string.h>
  8. #define random          rand
  9. #define bcopy(a, b, n)      _fmemmove(b, a, n)
  10. #define M_PI 3.14159265358979323846
  11.  
  12. #define MAXWINDOW        1000    /* Max analysis window length */
  13. #define FS                8000.0    /* Sampling rate */
  14.  
  15. #define DOWN            5        /* Decimation for pitch analyzer */
  16. #define PITCHORDER        4        /* Model order for pitch analyzer */
  17. #define FC                600.0    /* Pitch analyzer filter cutoff */
  18. #define MINPIT            50.0    /* Minimum pitch */
  19. #define MAXPIT            300.0    /* Maximum pitch */
  20.  
  21. #define MINPER            (int)(FS/(DOWN*MAXPIT)+.5)        /* Minimum period */
  22. #define MAXPER            (int)(FS/(DOWN*MINPIT)+.5)        /* Maximum period */
  23.  
  24. #define WSCALE            1.5863    /* Energy loss due to windowing */
  25.  
  26. static int framelen, buflen;
  27.  
  28. static float s[MAXWINDOW], y[MAXWINDOW], h[MAXWINDOW];
  29. static float fa[6], u, u1, yp1, yp2;
  30.  
  31. static int pitchctr;
  32. static float Oldper, OldG, Oldk[LPC_FILTORDER+1];
  33. static float b[LPC_FILTORDER+1], bp[LPC_FILTORDER+1], f[LPC_FILTORDER+1];
  34.  
  35. #define lin_to_ulaw(x) audio_s2u(x)
  36. #define ulaw_to_lin(x) audio_u2s(x)
  37.  
  38. static void
  39. auto_correl(float w[], int n, int p, float r[])
  40. {
  41.   int i, k, nk;
  42.  
  43.   for (k=0; k <= p; k++) {
  44.     nk = n-k;
  45.     r[k] = 0;
  46.     for (i=0; i < nk; i++) r[k] += w[i] * w[i+k];
  47.   }
  48. }
  49.  
  50. static void
  51. durbin(float r[], int p, float k[], float *g)
  52. {
  53.   int i, j;
  54.   float a[LPC_FILTORDER+1], at[LPC_FILTORDER+1], e;
  55.   
  56.   for (i=0; i <= p; i++) a[i] = at[i] = 0;
  57.     
  58.   e = r[0];
  59.   for (i=1; i <= p; i++) {
  60.     k[i] = -r[i];
  61.     for (j=1; j < i; j++) {
  62.       at[j] = a[j];
  63.       k[i] -= a[j] * r[i-j];
  64.     }
  65.     k[i] /= e;
  66.     a[i] = k[i];
  67.     for (j=1; j < i; j++) a[j] = at[j] + k[i] * at[i-j];
  68.     e *= 1.0 - k[i]*k[i];
  69.   }
  70.  
  71.   *g = sqrt(e);
  72. }
  73.  
  74. static void inverse_filter(float w[], float k[])
  75. {
  76.   int i, j;
  77.   float b[PITCHORDER+1], bp[PITCHORDER+1], f[PITCHORDER+1];
  78.   
  79.   for (i = 0; i <= PITCHORDER; i++) b[i] = f[i] = bp[i] = 0.0;
  80.     
  81.   for (i = 0; i < buflen/DOWN; i++) {
  82.     f[0] = b[0] = w[i];
  83.     for (j = 1; j <= PITCHORDER; j++) {
  84.       f[j] = f[j-1] + k[j] * bp[j-1];
  85.       b[j] = k[j] * f[j-1] + bp[j-1];
  86.       bp[j-1] = b[j-1];
  87.     }
  88.     w[i] = f[PITCHORDER];
  89.   }
  90. }
  91.  
  92. static void calc_pitch(float w[], float *per)
  93. {
  94.   int i, j, rpos;
  95.   float d[MAXWINDOW/DOWN], k[PITCHORDER+1], r[MAXPER+1], g, rmax;
  96.   float rval, rm, rp;
  97.   float a, b, c, x, y;
  98.   static int vuv=0;
  99.  
  100.   for (i=0, j=0; i < buflen; i+=DOWN) d[j++] = w[i];
  101.   auto_correl(d, buflen/DOWN, PITCHORDER, r);
  102.   durbin(r, PITCHORDER, k, &g);
  103.   inverse_filter(d, k);
  104.   auto_correl(d, buflen/DOWN, MAXPER+1, r);
  105.   rpos = 0;
  106.   rmax = 0.0;
  107.   for (i = MINPER; i <= MAXPER; i++) {
  108.     if (r[i] > rmax) {
  109.       rmax = r[i];
  110.       rpos = i;
  111.     }
  112.   }
  113.   
  114.   rm = r[rpos-1];
  115.   rp = r[rpos+1];
  116.   rval = rmax / r[0];
  117.  
  118.   a = 0.5 * rm - rmax + 0.5 * rp;
  119.   b = -0.5*rm*(2.0*rpos+1.0) + 2.0*rpos*rmax + 0.5*rp*(1.0-2.0*rpos);
  120.   c = 0.5*rm*(rpos*rpos+rpos) + rmax*(1.0-rpos*rpos) + 0.5*rp*(rpos*rpos-rpos);
  121.  
  122.   x = -b / (2.0 * a);
  123.   y = a*x*x + b*x + c;
  124.   x *= DOWN;
  125.  
  126.   rmax = y;
  127.   rval = rmax / r[0];
  128.   if (rval >= 0.4 || (vuv == 3 && rval >= 0.3)) {
  129.     *per = x;
  130.     vuv = (vuv & 1) * 2 + 1;
  131.   } else {
  132.     *per = 0.0;
  133.     vuv = (vuv & 1) * 2;
  134.   }
  135. }
  136.  
  137. int
  138. lpc_init(int len)
  139. {
  140.   int i;
  141.   float r, v, w, wcT;
  142.  
  143.   framelen = len;
  144.   buflen = len*3/2;
  145.   if (buflen > MAXWINDOW) return -1;
  146.  
  147.   for (i = 0; i < buflen; i++) {
  148.     s[i] = 0.0;
  149.     h[i] = WSCALE*(0.54 - 0.46 * cos(2 * M_PI * i / (buflen-1.0)));
  150.   }
  151.  
  152.   wcT = 2 * M_PI * FC / FS;
  153.   r = 0.36891079 * wcT;
  154.   v = 0.18445539 * wcT;
  155.   w = 0.92307712 * wcT;
  156.   fa[1] = -exp(-r);
  157.   fa[2] = 1.0 + fa[1];
  158.   fa[3] = -2.0 * exp(-v) * cos(w);
  159.   fa[4] = exp(-2.0 * v);
  160.   fa[5] = 1.0 + fa[3] + fa[4];
  161.  
  162.   u1 = 0.0;
  163.   yp1 = 0.0;
  164.   yp2 = 0.0;
  165.  
  166.   Oldper = 0.0;
  167.   OldG = 0.0;
  168.   for (i=1; i <= LPC_FILTORDER; i++) Oldk[i] = 0.0;
  169.   for (i=0; i <= LPC_FILTORDER; i++) b[i] = bp[i] = f[i] = 0.0;
  170.   pitchctr = 0;
  171.  
  172.   return 0;
  173. }
  174.  
  175. void
  176. lpc_analyze(unsigned char *buf, lpcparams_t *params)
  177. {
  178.   int i, j;
  179.   float w[MAXWINDOW], r[LPC_FILTORDER+1];
  180.   float per, G, k[LPC_FILTORDER+1];
  181.  
  182.   for (i=0, j=buflen-framelen; i < framelen; i++, j++) {
  183.     s[j] = ulaw_to_lin(buf[i])/32768.0;
  184.     u = fa[2] * s[j] - fa[1] * u1;
  185.     y[j] = fa[5] * u1 - fa[3] * yp1 - fa[4] * yp2;
  186.     u1 = u;
  187.     yp2 = yp1;
  188.     yp1 = y[j];
  189.   }
  190.  
  191.   calc_pitch(y, &per);
  192.  
  193.   for (i=0; i < buflen; i++) w[i] = s[i] * h[i];
  194.   auto_correl(w, buflen, LPC_FILTORDER, r);
  195.   durbin(r, LPC_FILTORDER, k, &G);
  196.  
  197.   params->period = (unsigned short) (per * (1<<8));
  198.   params->gain = (unsigned char) (G * (1<<8));
  199.   for (i=0; i < LPC_FILTORDER; i++) params->k[i] = (char) (k[i+1] * (1<<7));
  200.  
  201.   bcopy(s+framelen, s, (buflen-framelen)*sizeof(s[0]));
  202.   bcopy(y+framelen, y, (buflen-framelen)*sizeof(y[0]));
  203. }
  204.  
  205. int
  206. lpc_synthesize(lpcparams_t *params, double speed, unsigned char *buf)
  207. {
  208.   int i, j, flen=(int) (framelen/speed);
  209.   float per, G, k[LPC_FILTORDER+1];
  210.   float u, NewG, Ginc, Newper, perinc;
  211.   float Newk[LPC_FILTORDER+1], kinc[LPC_FILTORDER+1];
  212.  
  213.   per = (float)params->period / (1<<8);
  214.   G = (float)params->gain / (1<<8);
  215.   k[0] = 0.0;
  216.   for (i=0; i < LPC_FILTORDER; i++) k[i+1] = (float)params->k[i] / (1<<7);
  217.  
  218.   if (per == 0.0) {
  219.     G /= sqrt(buflen/3.0);
  220.   } else {
  221.     i = (int) (buflen / per);
  222.     if (i == 0) i = 1;
  223.     G /= sqrt((float)i);
  224.   }
  225.  
  226.   Newper = Oldper;
  227.   NewG = OldG;
  228.   for (i=1; i <= LPC_FILTORDER; i++) Newk[i] = Oldk[i];
  229.     
  230.   if (Oldper != 0 && per != 0) {
  231.     perinc = (per-Oldper) / flen;
  232.     Ginc = (G-OldG) / flen;
  233.     for (i=1; i <= LPC_FILTORDER; i++) kinc[i] = (k[i]-Oldk[i]) / flen;
  234.   } else {
  235.     perinc = 0.0;
  236.     Ginc = 0.0;
  237.     for (i=1; i <= LPC_FILTORDER; i++) kinc[i] = 0.0;
  238.   }
  239.     
  240.   if (Newper == 0) pitchctr = 0;
  241.     
  242.   for (i=0; i < flen; i++) {
  243.     if (Newper == 0) {
  244.       u = (random()/32767.0) * NewG;
  245.     } else {
  246.       if (pitchctr == 0) {
  247.         u = NewG;
  248.         pitchctr = (int) Newper;
  249.       } else {
  250.         u = 0.0;
  251.         pitchctr--;
  252.       }
  253.     }
  254.       
  255.     f[LPC_FILTORDER] = u;
  256.     for (j=LPC_FILTORDER; j >= 1; j--) {
  257.       f[j-1] = f[j] - Newk[j] * bp[j-1];
  258.       b[j] = Newk[j] * f[j-1] + bp[j-1];
  259.       bp[j] = b[j];
  260.     }
  261.     b[0] = bp[0] = f[0];
  262.     
  263.     buf[i] = lin_to_ulaw((int) (b[0] * 32768.0));
  264.     
  265.     Newper += perinc;
  266.     NewG += Ginc;
  267.     for (j=1; j <= LPC_FILTORDER; j++) Newk[j] += kinc[j];
  268.   }
  269.   
  270.   Oldper = per;
  271.   OldG = G;
  272.   for (i=1; i <= LPC_FILTORDER; i++) Oldk[i] = k[i];
  273.  
  274.   return flen;
  275. }
  276.